Finding the ciliary beating pattern with optimal efficiency 
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We introduce a measure for energetic efficiency of biological cilia acting individually or collectively 
and numerically determine the optimal beating patterns according to this criterion. Maximizing the 
efficiency of a single cilium leads to curly, often symmetric and somewhat counterintuitive patterns. 
But when looking at a densely ciliated surface, the optimal patterns become remarkably similar to 
what is observed in microorganisms like Paramecium. The optimal beating pattern then consists 
of a fast effective stroke and a slow sweeping recovery stroke. Metachronal coordination is essential 
for efficient pumping and the highest efficiency is achieved with antiplectic waves. Efficiency also 
increases with an increasing density of cilia up to the point where crowding becomes a problem. We 
finally relate the pumping efficiency of cilia to the swimming efficiency of a spherical microorganism 
and show that the experimentally estimated efficiency of Paramecium is surprisingly close to the 
theoretically possible optimum. 
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Many biological systems have evolved to work with a 
very high energetic efficiency. For example, muscle can 
convert the free energy of ATP hydrolysis to mechani- 
cal work with more than 50% efficiency ^Lj, the Fl-FO 
ATP synthase converts electrochemical energy of pro- 
tons to chemical energy stored in ATP molecules with 
even higher efficiency [2 , etc. At first glance, the beat- 
ing of cilia and fiagella does not fall into the category 
of processes with such a high efficiency. Cilia are hair- 
like protrusions that beat in an asymmetric fashion in 
order to pump the fiuid in the direction of their effec- 
tive stroke [3 . They propel certain protozoa, such as 
Paramecium^ and also fulfill a number of functions in 
mammals, including mucous clearance from airways, L- 
R asymmetry determination and transport of an egg cell 
in Fallopian tubes. Lighthill [4 defines the efficiency of 
a swimming microorganism as the power that would be 
needed to drag an object of the same size with the same 
speed through viscous fluid, divided by the actually dis- 
sipated power. Although the efficiency defined in this 
way could theoretically even exceed 100% [5] , the actual 
swimming efficiencies are of the order of 1% [6", 'T. In his 
legendary paper on life at low Reynolds number [8 Pur- 
cell stated that swimming microorganisms have a poor 
efficiency, but that the energy expenditure for swimming 
is so small that it is of no relevance for them (he uses 
the analogy of "driving a Datsun [a fuel-efficient car of 
the period] in Saudi Arabia"). Nevertheless, later studies 
show that swimming efficiency is important in microor- 
ganisms. In Paramecium^ more than half of the total 
energy consumption is needed for ciliary propulsion [9 . 

When applied to ciliary propulsion, Lighthill's effi- 
ciency [4 has some drawbacks. For one, it is not a direct 
criterion for the hydrodynamic efficiency of cilia as it also 
depends on the size and shape of the whole swimmer. Be- 
sides that it is, naturally, only applicable for swimmers 
and not for other systems involving ciliary fiuid transport 
with a variety of functions, like L-R asymmetry determi- 



nation [W. We therefore propose a different criterion 
for efficiency at the level of a single cilium or a carpet of 
cilia. A first thought might be to define it as the volume 
flow rate of the transported fluid, divided by the dissi- 
pated power. However, as the flow rate scales linearly 
with the velocity, but the dissipation quadratically, this 
criterion would yield the highest efficiency for infinitesi- 
mally slow cilia, just like optimizing the fuel consumption 
of a road vehicle alone might lead to fitting it with an 
infinitesimally weak engine. Instead, like engineers try 
to optimize the fuel consumption at a given speed, the 
well-posed question is which beating pattern of a cilium 
will achieve a certain fiow rate with the smallest possible 
dissipation. 

The problem of finding the optimal strokes of hypo- 
thetical microswimmers has drawn a lot of attention in 
recent years. Problems that have been solved include 
the optimal stroke pattern of Purcell's three link swim- 
mer [IT], an ideal elastic fiagellum [12 , a shape-changing 
body [13], a two- and a three-sphere swimmer [14 and a 
spherical squirmer [5] . Most recently, Tam and Hosoi op- 
timized the stroke patterns of Chlamydomonas^ fiagella 
[15 . But all these studies are still far from the complex- 
ity of a ciliary beat with an arbitrary three-dimensional 
shape, let alone from an infinite field of interacting cilia. 
In addition, they were all performed for the swimming 
efficiency of the whole microorganism, while our goal is 
to optimize the pumping efficiency at the level of a single 
cilium, which can be applicable to a much greater variety 
of ciliary systems. 

So we propose a cilium embedded in an infinite plane 
(at z = 0) and pumping fiuid in the direction of the 
positive X-axis. We define the volume fiow rate Q as the 
average fiux through a half-plane perpendicular to the 
direction of pumping [16 . With P we denote the average 
power with which the cilium acts on the fiuid, which is 
identical to the total dissipated power in the fiuid filled 
half-space. We then define the efficiency in a way that is 



independent of the beating frequency cj as 



(1) 



As we show in Appendix 1, minimizing the dissipated 
power P for a constant volume flow rate Q is equivalent 
to maximizing e at a constant frequency. A similar argu- 
ment for swimming efficiency has already been brought 
forward by Avron et al. [13 . 

Furthermore, a general consequence of low Reynolds 
number hydrodynamics is that the volume flow only de- 
pends on the shape of the stroke and on the frequency, 
but not on the actual time dependence of the motion 
within a cycle. This is the basis of Purcell's scallop the- 
orem [51. As a consequence, the optimum stroke always 
has a dissipation rate constant in time. We show this in 
Appendix 2. 

We can make the efficiency e completely dimensionless 
if we factor out the effects of the ciliary length L, the 
beating frequency uo and the fluid viscosity t]. The veloc- 
ity with which a point on the cilium moves scales with 
ujL and the linear force density (force per unit length) 
with rjUjL. The total dissipated power P, obtained by 
integration of the product of the velocity and linear force 
density over the length, then scales with rjuj^L^. The 
volume flow rate Q scales with ujL^. Finally, the effi- 
ciency e scales with L^ /rj. The dimensionless efficiency 
can therefore be defined as 



t' = L- 



r]e 
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When optimizing the efficiency of ciliary carpets, we 
have to use the measures of volume flow and dissipation 
per unit area, rather than per cilium. We introduce the 
surface density of cilia p, which is l/cP on a square lattice. 
In the following we show that the volume flow generated 
per unit area, pQ, is also equivalent to the flow veloc- 
ity above the ciliary layer. The fluid velocity above an 
infinite ciliated surface namely becomes homogeneous at 
a distance sufficiently larger than the ciliary length and 
metachronal wavelength. The far field of the flow induced 
by a single cilium located at the origin and pumping fluid 
in direction of the x axis has the form p/Tj 



v{x,y,z) = A—Cr 



(3) 



with an arbitrary amplitude A. For this field the volume 
flow rate is 
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dzv:^{x,y,z) = -A (4) 



and the velocity above an infinite field of such cilia is 

/CX) /»00 r) 

dx / dypvx{-x, -y, z) = -—pA = irpQ , 
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which is independent of z. In this regime, one can sim- 
plify the description of cilia by replacing them with a 
surface slip term with velocity Vc [H] . 

We now define the collective efficiency as Cc = 
{pQY / {pP) and in dimensionless form as 
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e^ is a function of the beat shape, the dimensionless den- 
sity pL? and the metachronal coordination, which will be 
explained later. Additionally, for a single cilium or for 
collective cilia the efficiency also depends on the dimen- 
sionless radius of the cilium, a/L, but this dependence is 
rather weak, of logarithmic order. 

At this point we note that our definition of efficiency is 
different from that used by Gueron and Levit-Gurevich 
[19]. They define efficiency as volume flux through a 
specifically chosen rectangle above the group of cilia di- 
vided by the dissipated power. While this measure is 
useful for studying the effect of coupling and metachronal 
coordination (they show that the collective efficiency of 
a group of cilia increases with its size), it lacks the scale 
invariance discussed above. Ganger et al. [20] studied a 
model for individual and collective magnetically driven 
artificial cilia. Rather than introducing a single measure 
for the efficiency, they studied the pumping performance 
(which is the more relevant quantity in artificial systems) 
and dissipation separately. They showed that the pump- 
ing performance per cilium can be improved with the 
proper choice of the metachronal wave vector, while the 
dissipation per cilium remains largely constant. Both 
studies were limited to two-dimensional geometry (pla- 
nar cilia arranged in a linear row) and neither of them 
uses a scale-invariant efficiency criterion proposed here. 
On the other hand, Lighthill's criterion for swimming or- 
ganisms shares the same scaling properties as ours (it 
scales with the square of the swimming velocity, divided 
by dissipation), but differs in definition because it mea- 
sures the swimming and not the pumping efficiency. At 
the end we will show how the two measures are related 
to each other for a spherical swimmer. 

Our goal is to find the beating patterns that have the 
highest possible efficiency for a single cilium, as well as 
the beating pattern, combined with the density and the 
wave vector that give the highest efficiency of a ciliated 
surface. 



THE MODEL 

We describe the cilium as a chain of A^ touching beads 
with radii a. The first bead of a cilium has the center 
position xi = (0, 0, a), and each next bead in the chain is 
located at x^+i = x^ + 2a(sin di cos (\)i , sin di sin (\)i , cos di). 
The maximum curvature of the cilium is limited by the 



condition 



(Xi+i -Xi)(Xi -Xi_i) > {2af COS Prr 



(7) 



Naturally, beads cannot overlap with the surface {zi > a) 
or with each other |xi — Xj| > 2a. 

We describe the hydrodynamics using the mobility ma- 
trix formalism. If the force acting on bead i is F^, the 
resulting velocities are 



dt 



tX,; 



E^^'^-F^- 



(8) 



In this formalism, each element Mij is itself a 3 x 3 ma- 
trix, corresponding to 3 spatial dimensions. In general, 
the above equation should also include angular veloci- 
ties and torques, but they are negligible for small beads 
when the surface speeds due to rotational motion are 
much smaller than those due to translational motion. 
The mobility matrix is symmetric and positive-definite 
[2r. Therefore, one can always invert it to determine the 
friction matrix T = M~^, which determines the forces on 
particles moving with known velocities 



Er.. 



(9) 



If the particles were at large distances relative to their 
sizes, the elements of the mobility matrix would be deter- 
mined by Blake's tensor [22 , which represents the Green 
function of the Stokes fiow in the presence of a no-slip 
boundary. In our case the condition of large interparti- 
cle distances is not fulfilled and we use the next higher 
approximation, which is the Rotne-Prager tensor in the 
presence of a boundary, as described in a previous paper 
[23. 

The volume fiow rate in x direction, averaged over one 
beat period T, depends on x-components of forces acting 
on particles and their heights z above the boundary [16 : 
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Zi{t)F^^i{t)dt . 



(10) 



The dissipation rate is simply the total power needed to 
move the beads against viscous drag. 



E--F' 



(11) 



We numerically maximized the quantity Q^ / P for a 
set of angles P^nax and different numbers of beads. We 
used the sequential quadratic programming algorithm 
(SQP) from NAG numerical libraries (Numerical Algo- 
rithms Group). The full details of the numerical proce- 
dure are given in Appendix 3. 

To study the collective efficiency and metachronal co- 
ordination, we studied an array of Na x Na cilia (unit 
cell) on a square lattice with a lattice constant d. We 



introduced periodic boundary conditions by adding hy- 
drodynamic interactions between particles and the rep- 
resentations beyond lattice boundaries. So if a certain 
element in the mobility matrix describing interaction be- 
tween particles at x^ and Xj is Mi^j(xi,Xj), we replace it 



byM;,,.(x„x,) = E^ 



^Mi^j(yii,yij -^pAcx -^qAcy). 



Here A = Nad denotes the size of the unit cell. For 
the sake of numerical efficiency, we used the full Rotne- 
Prager form for the first o instances (p, q = — o, . . . o) and 
approximated the interaction with its long range limit, 
independent of the actual particle positions, for the rest 
(SI). 

We expect the optimal solution to have the form of 
metachronal waves with a wave vector k = {kx-,ky) = 
{2t: /A){i<ix^ i^y). In order to satisfy the periodic boundary 
conditions, f<ix and tvy have to be integer numbers, e.g., 
between and Nn — 1. 



RESULTS 

Single-particle model 

We first start with some simple models that are not 
necessarily feasible in practice, but allow important in- 
sight into how the optimum is achieved. We will follow 
the spirit of the model used to study the synchronization 
of cilia [17 , where we replace the cilium by a small spher- 
ical particle. There are many swimmer models building 
on similar assumptions, for example the three-sphere- 
swimmer [24], and they all have in common that they 
assume the connections between spheres to be very thin 
and neglect any hydrodynamic forces acting on them. 

So the first hypothetical model we study is a single 
sphere of radius a that can move along an arbitrary path 
x(a;t) in the half space above the boundary, but in order 
to mimic the tip of a cilium it has to stay within the 
distance L of the origin, |x| < L — a. In order to sim- 
plify the calculation we also assume that the sphere is 
small, a <^ L. In this limit, we can neglect the effect of 
the boundary on the hydrodynamic drag, which is then 
always 7 = Qirrja. The dissipated power is then simply 
P = 7x^. Because it has to be constant in time, we can 
also write it as 



P = -ff/T^ 



(12) 



with £ denoting the total distance traveled within one 
cycle and T its period. The average volume ffow follows 
from Eq. ( 10 ) as 

z{t)-fx{t)dt = ^ i zdx= — ^ (13) 
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where S is the area of the particle trajectory, projected 
onto the x — z plane. The resulting efficiency is (IT]) 



(14) 
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FIG. 1: Optimal trajectories of the 1-particle model. A) Idealized case of a small particle restricted to \x\ < L. The solution 
consists of piecewise circular arcs, determined by geometric parameters a and r/L. B) Numerical solutions for finite-sized 
particles, plotted for different ratios a/L. C) Optimal path for a particle at a constant distance from the origin, |x| = L — a, 
with a = 0.1 L. The transparent hemisphere symbolizes the surface on which the particle can move. D) Dimensionless efficiency 
e as a function of the dimensionless particle radius a/L. The black line shows the model with variable distance and the red 
line with a fixed distance from the origin. The dashed line shows the limit of small radii (a <^ L), e — 0.192 a/L. 



To find the optimal path, we thus have to maximize the 
area-to-circumference ratio of the path, while fulfilling 
the constraints z > and |x| < L. Obviously, there is no 
benefit in going out of the x — z plane, but there is cost 
associated with it. Therefore, the optimum trajectories 
will be planar. As any curve that minimizes its circumfer- 
ence at a fixed surface area, the unconstrained segments 
of the trajectory have to be circle arcs. The curve has the 
shape shown in Fig. [T]A. A numerical solution shows that 
the area-to-circumference ratio is maximal if the angle a 
defined in Fig.[TJ\ has the value a = 0.483. The resulting 
maximal efficiency in the limit a <C L is e = 0.192 L'^a/rj, 
or, in dimensionless form, e' = 0.192 a/L. 

Solutions for finite values of a/L are shown in Fig. [lb 
and their efficiencies in Fig.[T)l). The highest possible nu- 
merical efficiency of this model is e' = 0.0087, which is 
achieved at a/L = 0.13. 

Another version of the single-particle model is one in 
which the particle has to maintain a fixed distance (L — a) 
from the origin, while it is free to move along the surface 
of a sphere (Fig. [ip). This is an additional constraint 
and can therefore only reduce the achievable efficiency. 
As shown by the red line in Fig.[T)l), the efficiency indeed 
lies somewhat below that of the model with a variable 
distance and reaches a maximum value of e' = 0.0065. 



N particles, stiff cilium 

The next minimalistic model we will study is a stiff 
cilium: a straight chain of N beads with radius a and a 
total length of L = 2Na that can rotate freely around 



the center of the first bead. The problem is related to 
artificial cilia driven by a magnetic field [20l |23l |25] in 
which the orientation of the cilium largely (although not 
completely) follows the direction of the magnetic field. A 
related optimization has been performed by Smith et al. 
[16], but with two important differences. First, Smith et 
al. optimize the volume flow alone and not the efficiency. 
Their optimal stroke therefore touches the surface during 
the recovery stroke, while ours has to keep some distance 
in order to limit the dissipation. Second, they restrict 
themselves to cilia beating along tilted cones, whereas 
we allow any arbitrary pattern. 

The motion of a stiff cilium on its optimal trajectory is 
shown in Figure [2}^. The path of its tip closely resembles 
that of a single sphere at a fixed radius. The resulting 
dimensionless efficiency for A^ = 20 beads is e' ~ 0.00535. 



N particles, flexible cilium 

As the next level of complexity, and at the same time 
the first realistic description of biological cilia, we now 
study a flexible cilium consisting of A^ spherical particles 
(we use TV = 10 and N = 20). The bending angle per 
particle is restricted to /3max Q. Such a constraint is 
necessary for two reasons. For one, the curvature of a 
biological cilium is restricted by the bending rigidity of 
the axoneme. In addition, our A^-particle model verita- 
bly represents a continuous cilium only if the curvature 
radius is sufficiently larger than the size of a sphere. Ex- 
amples of beating patterns obtained by numerical opti- 
mization are shown in Figure [26,0. Figure [2t) shows the 
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dimensionless efficiency ^ as a function of /3max- 

It is instructive to look at fundamental symmetries of 
the problem at this point. First, as any of the problems 
studied here, it is symmetric upon reflection y -^ —y. For 
every clockwise beat, there is an equivalent counterclock- 
wise beat with the same efficiency. All cycles discussed 
here spontaneously break the y symmetry. In addition, 
the equations are invariant upon reflection x -^ —x with 
simultaneous time reversal, t -^ —t. This symmetry can 
be broken or not at the efficiency maximum. Interest- 
ingly, this depends on the allowed bending between adja- 
cent elements /3max- For example, the solution shown in 
Fig. |2b is xt-symmetric, while the one in Fig. [2p is not. 
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FIG. 2: Optimal beating patterns of a cilium consisting of 
N — 20 particles with different allowed bending angles /3max- 
The gray surface shows the projection of the tip trajectory on 
the X — y plane. A) A stiff cilium, /3max = 0. B) A flexible 
cilium, /3max = 20°. The optimal stroke is symmetric in x 
direction. C) Flexible cilium, /3max = 30°. The symmetry 
in X direction is broken. D) Dimensionless efficiency e as a 
function of /3max. The black line shows the optimal symmet- 
ric solution and the red line the asymmetric solution in cases 
where it is more efficient. The dashed line shows the maxi- 
mum efficiency for N = 10 and double /3max (corresponding 
to the same curvature). 



Multiple cilia and metachronal Avaves 

We solve the optimization problem of Na x A^^ cilia 
{Na = 12) with periodic boundary conditions by impos- 
ing a wave vector ikx^ky)^ finding the optimal solution 
for that vector and repeating the procedure for Na x Na 
wave vectors. Examples of optimal solutions for 6 differ- 
ent wave vectors are shown in Fig. [3]A.-F. The efficiency 
e'^ as a function of the wave vector is shown in Fig. [3J3. 
Note that all solutions determined in this section are for 
counterclockwise beating (viewed from above) . For clock- 
wise strokes the y component of the wave vector would 
have the opposite sign. Optimal solutions for four dif- 
ferent values of the interciliary distance d are shown in 
Figure |4|^-D and the optimal efficiency as a function of 
d in Figure [4^. 

Figure [3p shows that the efficiency depends more 
strongly on the longitudinal ikx) component of the wave 
vector than on the lateral {ky). This could partly be 
due to the nature of the hydrodynamic interaction, which 
is stronger in longitudinal direction, and partly because 
the cilia exert larger motion in longitudinal direction and 
therefore come closer to their neighbors along the x axis. 
Antiplectic metachronal waves (waves propagate in the 
opposite direction from the fluid pumping) generally have 
a higher efficiency than symplectic, but the fine structure 
is much more complex. For low ciliary densities, the op- 
timal solution is found for waves propagating along the x 
axis. However, for higher densities solutions with a pos- 
itive kx are more efficient, which means that the waves 
are dexio-antiplectic. Efficiency also grows with increas- 
ing density. But when the interciliary distance reaches 
d = 0.25 L crowding becomes a problem and the effi- 
ciency drops again. At even higher densities the solu- 
tion becomes increasingly difficult to find because of the 
complicated topology of densely packed cilia. Another 
problem is that the wavelength of the optimal solution, 
relative to the lattice constant, becomes increasingly long 
at high densities, which would require a unit cell larger 
than 12 x 12 used in our calculations. 






- ( — ( — ( — I — I — ( — r^ 



71^ 



rt^ 



rx- -^ rx^ 



cx^ 



"rx 



rx. - rx- 



rx^ 
rx. 



.— CX_ - 


^— CX_ - 


--ct 


^Ct 


^CX- ~- 


.— CX_ ~- 


^cx . 


^cx . 


^ct 


^CX_ -^ 


_cx. - 


_cx. _ 



o o o o 
r c r c 
r~ c~ r~ c- 



o o o o o 
r- r- r- r- r 



Cb 

- c 

- c 

- c 

- c 
. c 

- c 

- c 



c- 
c- 



f: 



r~ r- r- r- 
( t ( I 
o o o o 



r- o- ^ r~ r- c 

([III 

o o o -o -o 




0.007 



0.002 



FIG. 3: Optimal solutions at fixed wave vectors for interciliary distance d — 1.^ x L^ N — 20, /3max = 15° and Na — 12. A-F) 
Optimal solutions for wave vectors (fe,A:^) = (0,0) (A), (-7r/(2d), 0) (B), (57r/(6d),0) (C), (-27r/(3d), vr/d) (D), (0, -7r/(3(i)) 
(E) and (0,7r/(3(i)) (F). The blue arrow (x axis) denotes the direction of pumping and the red arrow the wavelength and 
direction of metachronal wave propagation. G) Efficiency e^ (red color represents high efficiency) as a function of the wave 
vector (kx^ ky). The maximum efficiency is in this case achieved for k = (— 7r/(2(i), 0) and antiplectic waves are generally more 
efficient than symplectic. The synchronous solution (0, 0) represents the global minimum of efficiency. 



Swimming efficiency of a ciliated microorganism 

We can finally use these results to estimate the maxi- 
mum possible swimming efficiency of a ciliated microor- 
ganism. For the sake of simplicity, we assume that the 
swimmer has a spherical shape with radius R. Accord- 
ing to Lighthill's definition, the swimming efficiency is 
defined as 



CL 



GtttjRV^ 



where V is the velocity and Ptot the total dissipated 
power [4^. Assuming that the layer of cilia is thin in 
comparison with the size of the organism {L <C R), the 
swimming velocity V can be calculated as 



V 



= 4^ /"^'^ =i£ .(^)2.sin^ OdO . (16) 



Here v{6) is the propulsion velocity above the ciliated 
layer. The dissipation can be expressed as 
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v^ (i9) sin (9d(9. 



multiplier 
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5v{e) 5v{e) 



(18) 



which is fulfilled if the angular dependence has the form 
v{0) = vq sinO. This leads to 
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(19) 



(15) and 



87^RX 



(20) 



Together, these equations give Lighthill's efficiency 



2 ^ / 



(21) 



With a maximum e^ ~ 0.016 and a typical ratio L/R ^ 
0.1, we obtain e^ ^ 0.016. 

For comparison, an optimized envelope model yields 
an efficiency around 3 % with the same parameters (if we 
translate the ciliary length into maximal displacement of 
the envelope) [5], which shows that the latter is a rela- 
tively good approximation for cilia if they operate close 
to the optimal way. 



In the second equality we used the definition (|6|, as well 
as the relationship between v and pQ ([s]). In order to 
obtain the maximum F at a given Ptot, the two func- 
tional derivatives should be related through a Lagrange 



DISCUSSION 

We introduced a scale invariant efficiency measure for 
the fluid pumping by cilia and started with optimizing 
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FIG. 4: Optimal solutions for various interciliary distances: d = 1.0 L (A), d — 0.7L (B), d — O.^L (C) and d — 0.25 L (D). For 
reasons of clarity the front rows of cilia are omitted and instead of individual spheres used in the calculation a tube connecting 
them is shown. E) Highest efficiency e^ as a function of the interciliary distance d for cilia consisting oi N — 20 (circles) and 
A^ = 10 (squares) spheres. For N = 10, we set /3max = 30° in order to allow the same maximum curvature. 



three simple instructive systems: a free sphere, a sphere 
at a constant distance from the origin and a stiff cilium. 
Of those three, the free sphere can reach the highest ef- 
ficiency. But they are all topped by a flexible cilium. 
The thickness of a cilium has only a small effect on its 
efficiency, which strengthens our choice to describe the 
cilium as a chain of beads. Depending on the allowed 
bending, the flexible cilium can have different shapes of 
the optimal beating pattern. In most cases the cilium 
curls up during the recovery stroke, rather than sweep- 
ing along the surface. Such shapes appear "unnatural" if 
we compare them with those observed in microorganisms 

But the collective optimization of ciliary carpets leads 
to beating patterns that are strikingly similar to what 
is observed in many ciliated microorganisms. Unlike iso- 
lated cilia, they contain a recovery stroke during which 
they sweep along the surface. This is primarily due to the 
fact that beating patterns that are optimal for a single 



cilium (e.g., as shown in Fig. ^p) are not possible on a 
dense grid due to steric hindrance. The sweeping recov- 
ery stroke, on the other hand, allows dense stacking of 
cilia (best seen in Fig. [4t)) which further reduces drag as 
well as backward flow. The optimal effective stroke be- 
comes significantly faster than the recovery stroke. While 
a single cilium reaches its highest efficiency if the effec- 
tive stroke takes about 45% of the cycle, the optimum is 
around 20 — 25% for densely packed cilia. A similar ra- 
tio has been observed in Paramecium [28]. The distance 
between adjacent cilia in Paramecium is between 0.151/ 
and 0.25 L [29j, consistent with the predicted optimum 
around d = 0.25 L. It is interesting that the efficiency 
of any other wave vector is higher than the efficiency 
of the synchronous solution (wave vector 0). This is in 
agreement with some previous simpler, one dimensional 
models [20], but has not yet been shown on a 2-D lat- 
tice. We also find that antiplectic waves are generally 
more efficient than symplectic, although symplectic solu- 



tions with a relatively high efficiency exist, too. For high 
densities and cilia beating counterclockwise, the waves 
become almost dexioplectic (meaning that the effective 
stroke points to the right of the wave propagation) and 
the wavelength becomes similar to the cilium length L 
- both findings are in agreement with observations on 
Paramecium [3l |30]. For cilia beating clockwise, laeo- 
plectic waves would be more efficient, which is indeed 
observed [31 . Although the effect of thickness is small, 
it is interesting to note that thicker cilia have a slightly 
higher efficiency when isolated or at low surface densities, 
but are outperformed by thinner cilia at high densities. 

The total energetic efficiency of swimming in Parame- 
cium has been measured as 0.078% [9]. This figure in- 
cludes losses in metabolism and force generation - the 
hydrodynamic swimming efficiency alone has been esti- 
mated as 0.77%. This comes close (by a factor of 2) to 
our result for the maximally possible Lighthill efficiency 
of a spherical ciliated swimmer cl ~ 0.016. A biflagellate 
swimmer like Chlamydomonas has a lower theoretical ef- 
ficiency of 0.008 [15], but it is still within the same order 
of magnitude. 

Although efficiencies below 1 % seem low, we have 
shown that Paramecium still works remarkably close to 
the maximum efficiency that can be achieved with its 
length of cilia. While longer cilia might have a higher 
swimming efficiency, there are other considerations that 
are not included in this purely hydrodynamic study. For 
example, the bending moments and the power output per 
ciliary length can be limiting [32 . Thus, our study shows 
that at least for ciliates like Paramecium^ Purcell's view 
that efficiency is irrelevant for ciliary propulsion has to be 
revisited. Efficiency of swimming does matter for them, 
and in their own world they have well evolved to swim 
remarkably close to the optimal way. 

We have benefited from fruitful discussions with Frank 
Jiilicher. This work was supported by the Slovenian Of- 
fice of Science (Grants Pl-0099, Pl-0192, Jl-2209, Jl- 
2200 and Jl-0908). 



proves immediately that minimizing the dissipated power 
while keeping the volume fiow Q constant is equivalent 
to maximizing e at a constant uj. 



Appendix 2. Optimal solutions have a constant 
dissipation 

In the following we demonstrate that a trajectory with 
optimal efficiency always has a dissipation that is con- 
stant in time. The pumping performance only depends 
on the shape of the trajectory and the period T, but not 
on the velocity along the trajectory. We have to deter- 
mine the latter in a way that minimizes the average dis- 
sipation. We parameterize the trajectory as x((/p) where 
the phase (/:? is a function of time that fulfills (/p(0) = and 
(/p(T) = 27r. The dissipation at any moment is quadratic 
in velocity, P{t) = x-^7(x)x, and can be written as 

P{t) = {d^/diff-f{^){d^/dif)ip^ = p{ip)if^ . (S2) 

We obtain the average dissipation by integrating over one 
period, 

rT 1 rT 



P{t)dt 



1 
f 



p{(p)(p dt 



, M^ 



dip ■■ 



1 
T 



p{if,t')dif , 



(S3) 

with t' = dt/d(p and p{(p,t') = p{(p)/t'. Minimization of 
the average dissipation while keeping the period constant 
leads to the following Euler-Lagrange equation (note the 
swapped roles of t and (p) 

dp d dp dp dP 

m ~ dApdf ~ dAft^ ~ ~dAp ' 

This proves that the optimal solution has a dissipation 
P that is constant in time. 







(S4) 



Appendix 3. Numerical optimization procedure 



Appendix 1. Justification of the efficiency criterion 

In this section we prove that instead of minimizing the 
dissipated power P for a constant volume fiow rate Q 
we can maximize the numerical efficiency e = Q^ /P at a 
constant beating frequency uj. Let the volume flow rate 
Q[x(a;t)] and dissipation P[x(a;t)] be functionals of the 
trajectory shape x(cjt) and functions of cj. The efficiency 
e, defined as e[x(a;t)] = Q^/P^ is only a functional of 
x(cc;t), but independent ofuj. The relationship 



Q^ SP 
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Single cilium 

For a single cilium, modeled as a chain of N beads, 
we parameterize the stroke as a sequence of Ns equally 
spaced time steps with duration At = T/Ns- Let the 
vector x^(r) represent the position of the bead i at time 
step r. Because the trajectory is periodic in time, we 
have x^(r) = x^(r + Ns). For a given stroke the volume 
fiow rate is calculated according to Eq. [lO] as 

Q^^^Y^Y^ ^.(r + l)r,,,(r + l) + z,(r)r,,,(r) 



TTrjNs 



r=l i,jf = l 



,Xj(t+1)-Xj(t) 



M 



(S5) 



T] denotes the viscosity of the surrounding fluid and Zi 
the height of i-th bead above the surface. The dissipated 
power, defined in Eq. [TT] can be written as 

-, Ns N 

=ii,i=i (S6) 



r,,,(r + i) + r,,,(r) 



'Ns 



(x,(t + 1) - Xi(r)) (xj(r + 1) - Xj(r)) 
(Ai)2 

The dimensionless efficiency, which we want to maximize, 
follows as e' = L-^rjQ^/P with L = 2aN ^. The fric- 
tion matrix Tij{t) that appears in above expressions is 
obtained by numerical inversion of the mobility matrix 
Mij^ which is calculated using the Rotne-Prager approx- 
imation in the presence of a boundary, as described in 

We computed the optimal strokes using the NAG (The 
Numerical Algorithms Group Ltd., Oxford, UK) Fortran 
Library E04UGF routine, which maximizes an arbitrary 
smooth function subject to constraints using a sequen- 
tial quadratic programming method. At each time step 
we parameterize the bead positions x^ by a set of A^ — 1 
pairs of polar {0) and azimuthal ((/>) angles, as described 
in the main text (Section "The Model"). Therefore, the 
parameter space in which we search for the for the effi- 
ciency maximum is 2Ns{N — 1) dimensional. For most 
calculations presented here {N = 20, Ns = 84), this 
means 3192 dimensions. We still consider optimization 
in such high-dimensional space preferable to parameter- 
izing the beat shape in terms of Fourier modes (see, e.g. 
[15]), which would have less variables, but therefore a 
more complex landscape with more local maxima. 

The stroke shapes are subject to two types of con- 
straints: (i) beads are modeled as hard core objects that 
allow no overlapping of two beads or a bead with the 
surface, (ii) the maximum bending angle defined between 
three adjacent beads is limited to /3max 0. 

We initialized the optimization procedure with a sim- 
ple beating pattern (a tilted cone) that had the correct 
handedness (clockwise) and checked that the result was 
otherwise independent of the initial state. 



Carpets of cilia 

For an array of cilia the size of parameter space would 
grow with the square of the system size Na and the 
demand to calculate the hydrodynamic mobility matrix 
with its fourth power. The problem would become nu- 
merically unsolvable even for a relatively small field of 
cilia. However, in an infinite system (or a system with 
periodic boundary conditions) we expect the optimal so- 
lution to have the form of metachronal waves with an 
unknown wave vector. The optimization problem at a 
fixed wave vector then requires the same number of pa- 
rameters as the optimization of a single ciliary beat. The 
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Long-range approximation 



FIG. SI: The elements of the mobility matrix describe the 
response of the particles forming the central cilium (0, 0) to 
forces acting on all other particles. Cilia are positioned on 
a square lattice with lattice constant d, while the unit cell 
comprising Na x Na cilia has the edge length A. To facilitate 
the calculation we use the full Rotne-Prager form for (2o + 
1) X (2o+l) unit cells and use the long-range limit expressions 
for more distant cilia. 



globally optimal solution can be found by repeating the 
calculation for all wave vectors that fulfill the periodic 
boundary conditions. 

For a wave vector {kx^ky)^ the i-th sphere of the cilium 
at lattice position (a, /3) follows the path 

Xi,a,/3(cjt) = adex-\-f^dey-\-:sii{ujt — adkx — pdky) , (S7) 

where x^(co't) is the path of the cilium at lattice position 
(0,0), which we are optimizing. Similarly, the force act- 
ing on the same bead has the following time dependence 



Fi^a,i3{(^t) = Fi{ujt -adkx- pdky) . 



(S8) 



The equations of motion of z-th bead belonging to the 
cilium (0, 0) read 

N oo oo 

= 2_\ /_! Z_l Mij.oc^j3{ujt)¥j{ujt — adkx — pdky) 

jf = l Q! = — CX) /3 = — OO 

(S9) 

In this expression Mij.a,/3{^t) represents the element of 
the mobility matrix that links the response of i-th bead 
of the cilium (0, 0) to the force acting on j-th bead of the 
cilium (a, /3) at time t. 
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Periodic boundary conditions are fulfilled if kx and ky 
are multiples of 27r/A, with A = Nad denoting the size 
of the field (Fig. [sT]) . We can therefore use integer wave 



vectors Hix^y = kx^yA/{27T) instead. We can make use 
of the periodicity by carrying out the summation in the 
above equation 



d 
~df 



N 



-^iiujt) = ^^Mlj.^^^{ujt)Fj{ujt - adkx - pdky) 

(SIO) 
with 



m;, 



i,j;Q;,/3 — 2-^ Mij.a-\-pNa 



.^+qNa 



(Sll) 



The indices a and /3 run over one unit cell, centered 
around the origin. For odd A^o, this would be from 
-{Na - l)/2 to {Na - l)/2. For an even A^^, they run 
from —Na/2 to Na/2^ but the boundary terms are given 
half weight each. This distribution is necessary in order 
to preserve the symmetry of the mobility/friction matrix. 
For the sake of numerical efficiency, we use the full 
Rotne-Prager form for the first o instances (p, q = 
—o^...o) and approximate the interaction with 



({pA^adY {pA^ad){qA^ pd) 

[pA + ad) [qA + pd) [qA ^ pdf 



(S12) 

for larger indices. We used o = 1 for the optimization 
and checked the final efficiency with o = 2 (the correc- 
tion was not significant). The contribution of distant 
cilia towards M' can easily be calculated numerically to 
quadratic order in a and /3 and has the form 



27TT]A^ 



Ci 



C2 



Ni 



■C,^ 



3 7V2 



^4]v2 

^ 



C4^ 0\ 

d ^0/ 

(S13) 
Numerical values of the coefficients Ci . . . C4 are listed in 
table EH 



TABLE SI: Coefficients C1...C4 describing the contribu- 
tions of cilia beyond o periods, for different orders o. 



o Ci 



Co 



C3 



C4 



4.51681 14.0613 -2.60821 -12.8518 

1 1.80961 0.680011 0.182076 -0.210573 

2 1.1135 0.161671 0.0474549 -0.0445071 

3 0.801452 0.0608975 0.0181854 -0.0163511 

With discretized time steps the expression for velocity 



(SIO) becomes 



Xi(r) 



N 

E E K.;«,/3(^)F, (r - aNfH^x - fSNfKy) 



N Ns 



(S14) 



with M,^j.r,r' = Ea,^Kr.aA^)l^' = ^ ' ^^/^- ' 

PNfKy (mod Ns)] denoting the generalized mobility ma- 
trix that now couples the forces and velocities at different 
time steps. The brackets [] denote a function which is 1 if 
the condition is fulfilled and otherwise. We also set the 
number of time steps to be a multiple of the lattice size, 
Ns = NfNa. M is in total a '^NNs x 'iNNs dimensional 
matrix, which can be re-ordered into a block-diagonal 
form with Nf blocks for greater numerical efficiency. 

We calculate corresponding generalized friction matrix 
by inverting the generalized mobility matrix. 



M-^ . 



(S15) 



The friction matrix now allows us to calculate the forces 
if we know the velocities of all beads at all times. 



N Ns 



(S16) 



Within the finite difference approximation this leads to 
equations 



1 






r,r' = l i,j = l 
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(S17) 
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1 V^ V^ ri,j;r+l 

N^ ^ ^ 



,r' + l ~r i i,j;T,T' 



t,t' = 1 ij = l 

(x,(r' + l)-x,(r'))(x,(T' + l)-x,(r')) 



(S18) 



(Aty 



Now we can numerically optimize the quantity QjyP as 
a function of the angles that parameterize the beat shape 
(in total 2(A^ — \)Ns variables). We initialize the system 
in a way that only solutions with clockwise rotation (as 
seen from above) are considered. The constraints are 
similar to those on a single cilium, except that hard-core 
repulsion between neighboring cilia has to be taken into 
account, too. 

In conclusion, the combination of an efficient descrip- 
tion of hydrodynamics, the usage of periodic bound- 
ary conditions and the translational invariance of the 
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metachronal wave allowed us to reduce the optimization 
problem to one with the same number of variables as 
for a single cilium. The computational demand to calcu- 
late the efficiency for a certain set of variables scales as 
((2o+ l)A/'a)^, rather than the fourth power which would 
be the case if all cilia were treated independently. The 
optimization for one wave vector typically takes a few 
days on one core of the Intel Xeon E5520 CPU. To find 
the global maximum we have to repeat the calculation 
with all wave vectors, N'^ in total. 
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